UTHEP-284 
September 1994 



Critical exponents 
of a three dimensional 0(4) spin model 

K. Kanaya and S. Kaya 

Institute of Physics, University of Tsukuha, Ibaraki 305, Japan 

Abstract 

By Monte Carlo simulation we study the critical exponents governing the tran- 
sition of the three-dimensional classical 0(4) Heisenberg model, which is considered 
to be in the same universality class as the finite-temperature QCD with massless 
two flavors. We use the single cluster algorithm and the histogram reweighting 
technique to obtain observables at the critical temperature. After estimating an 
accurate value of the inverse critical temperature K^. = 0.9360(1), we make non- 
perturbative estimates for various critical exponents by finite-size scaling analysis. 
They are in excellent agreement with those obtained with the 4— e expansion method 
with errors reduced to about halves of them. 
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1 Introduction 



Finite temperature chiral phase transition of QCD is very important in the study of 
phase transitions in the early Universe and in the investigation of heavy ion colhsions 
at high energy. At present, this transition is studied mainly using the Monte Carlo 
method on the lattice. Pisarski and Wilczek ^ suggested that QCD with massless 
two flavors, which is considered to be an approximation of the real world, belongs to 
the same universality class as three-dimensional four-component Heisenberg models, 
if the finite temperature chiral transition of Nf = 2 QCD is second order. Then, 
the chiral transition of Nf = 2 QCD has the same critical exponents as the 3d 0(4) 
Heisenberg model. 

Simulations of lattice QCD for Nf = 2 suggest that the chiral transition is a 
second order transition for staggered fermions and for Wilson fermions Q. The 
study towards a precise measurement of the critical exponents of Nf = 2 QCD has 
just begun In the verification that the 0(4) Heisenberg model belongs to the 
same universality class, there is a problem that both Wilson fermions and staggered 
fermions on the lattice do not have the full chiral symmetry — which is expected to 
restore only in the continuum limit. Conversely, however, we could consider that, 
assuming the universality, the chiral symmetry is restored on the lattice sufficiently 
when the exponents agree with those of the 3d 0(4) Heisenberg model. 

Therefore an accurate calculation of the critical exponents of the 3d 0(4) Heisen- 
berg model is quite important. For this model the best estimation of critical expo- 
nents has been made with the 4 — e expansion method up to seven loops ||6|. 

In this work we simulate the 3d 0(4) Heisenberg model by the Monte Carlo 
method and make a non-perturbative estimation of several critical exponents. We 
use the single cluster Monte Carlo update algorithm which recently has been used 
for the simulation of spin systems: Wolff formulated this algorithm by modifying the 
multiple-cluster algorithm by Swendsen and Wang[^ and applied it to continuous 
spin models ^. Recent applications of the multiple and single cluster algorithms 
to two and three dimensional spin models have demonstrated their advantage in the 
computation time to the usual local update algorithms. Among global algorithms, 
the single-cluster algorithm is shown to be superior to the multiple-cluster algorithm 
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for three-dimensional spin models [10, Therefore, we apply the single-cluster 
algorithm in this study. 

In section |2|the model and the method of simulation are described. In section ^ 
we estimate the transition point from the crossing point of the Binder cumulant and 
compute the critical exponents at the transition point making use of the histogram 
reweighting technique. We also check the consistency of the results by independent 
measurements of the critical temperature and several exponents. We then compare 
our exponents with those of the 4 — e expansion method. Our conclusion is given 
in section ^. 



2 The model and the method 

The partition function Z and the energy E of the 3d 0(4) Heisenberg model are 
defined by 

z = n/['^^>™' 

X •' 

E = Y,{l-s{x)-s{x + i)}, (1) 
x,i 

where K is the inverse temperature and s{x) is a four-dimensional unit spin at the 
lattice site x. Vs are the unit steps in three coordinate directions. We use three- 
dimensional simple cubic lattices with the volume V = with L = 8, 10, 12, 14, 
16, 24, and 32, and employed periodic boundary conditions. 

We chose two simulation points for each L except for L = 10: One \s K = 0.935 
which is a rough estimate for the transition point by our preparatory simulation. 
Another simulation point is chosen for each L at the maximum point of the sus- 
ceptibility estimated by a preparatory simulation. Our simulation parameters are 
compiled in Table |^. We use the data at K = 0.935 for the calculation of the transi- 
tion point as well as the analyses of finite-size scaling with the histogram reweighting 
technique and use the data at the maximum of the susceptibility for a check of the 
consistency of our results. 

The magnetization and the energy are measured every 10 sweeps and stored on 
the disk. We define one sweep by one cluster update by the single-cluster algorithm 
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explained in the following section. Several million sweeps are performed for each 
simulation point. From the autocorrelation time we measured (see the following 
section) this corresponds to about one hundred thousand independent data for each 
point, as compiled in Table We estimate errors by the jackknife procedure. We 
study the bin-size dependence of erros and choose a suffiently large bin-size such 
that errors become stable. The resulting bin-sizes are consistent with the values 
of auto-correlation time estimated independently. All the jobs takes 23 hours with 
HITAC S820/80. 



We use the histogram reweighting method] 12 1 to calculate the observables in 
a region of K around the simulated point Kgimu- The region of K in which the 
histogram reweighting method is applicable can be determined by the magnitude 
of the shift of energy value: If the peak position of a reweighted energy distribu- 
tion, EpeakiK), locates away from the peak position of the original distribution, 
Epeak{K simu) 1 then the statistical errors for averages computed with the reweighted 
distribution become large correspondingly. Limited statistics near the tails of mea- 
sured histograms also leads to a danger of large under-estimation of the errors there. 
We study the effect of reweighting and observe that, with our statistics, many er- 
rors for the observables we study become rapidly large and the histogram becomes 
rapidly notched when K gets outside the region where the height of the original 
energy histogram at Ep(,ak{K) is larger than one third of the peak height. Although 
several computed errors, such as the error for the Binder cumulant discussed below, 
sometimes remain small even outside this range, we find that the result is not con- 
sistent with the result of a direct simulation there. We therefore limit ourselves to 
apply the histogram reweighting method only up to the point where the height of 
the original energy histogram at EpeakiK) decreases to a third of the peak height. 
Similar criterion is used also in Ref. |Tc 

2.1 Algorithm 

We use the single-cluster algorithm formulated by Wolff This is a global update 
algorithm whose advantage is that the autocorrelation time and the dynamical ex- 
ponent are both much smaller than those of the local update algorithm as discussed 
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below. 

The autocorrelation function A{k) is defined by 

p{k) = {OiO,+k) - {Oif, (2) 

with Oi being the i-th measurement of an observable O. The autocorrelation time 
r which is given by integrating the autocorrelation function 

oo 

r = Y.Aik) (3) 

fc=i 

diverges in the critical region as r oc where ^ is the correlation length. The 
exponent z is called as dynamical exponent. On finite lattices in the critical region, 
^ is replaced by the lattice length L: 

T(xL' . (4) 

This lattice size dependence of r is the origin of the "critical slowing down" which 
makes difficult to get a high effective statistics in the critical region on large lattices. 
We should use an algorithm with a small dynamical exponent. It is known that the 
local update algorithms such as the Metropolis algorithm have z ~ 2 independent 
of the model and the details of the update algorithm. For example, z = 1.94(6) is 
obtained for the 3d 0(3) Heisenberg model with Metropolis algorithm Use of 



a global update algorithm is required to get a smaller z. It is reported in Refs.|lC, 



13, ^ that z with the single-cluster algorithm for the 3d Ising model is about 0.2 
and that for the 3d 0(3) Heisenberg model is about 0. As presented in the next 
section, our result of z for the 3d 0(4) Heisenberg model is also consistent with 0. 
The single-cluster update for the 0(n) Heisenberg model is as follows |8|, [T6[ |: 

(1) A unit vector in the 0(n) space, r, is chosen with a random direction. 

(2) A starting site of a cluster, xq, is chosen at random and is included in the 
cluster. 

(3) For a link on the surface of the cluster, dc = {{x,y)\x & C , y ^ C)}, y is 
included in the cluster with the probability 

P{sx,sy) = 1 - exp[min{0, -2K{r- Sx){r ■ Sy)}] . (5) 
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(4) The process (3) is repeated until the growth of the cluster stops. 

(5) All spins in the cluster is flipped with regards to the surface perpendicular to 
f: 

s'x = R{r)sx = sx- 2{sx ■ r)f . (6) 

In order to test the efficiency of the algorithm and to test our program code 
for the single-cluster update, we simulate the 3d Ising and the 3d 0(3) Heisenberg 



model. Our results are completely consistent with Refs.[|Tn, 13], including the results 



for susceptibility, dynamical exponent, and critical exponents. 

3 Results 

3.1 Autocorrelation time and energy distribution 

Our results for the autocorrelation time are compiled in Table ||. The autocorrela- 
tion time Tra in terms of Metropolis unit stays almost constant or rather decreases 
with the increase of the lattice size. This implies that the dynamical exponent is 
or slightly smaller than 0. Similar result is obtained also for the 0(3) Heisenberg 
model [p^]. 

The measured energy distribution shown in Fig. |l] is a Gaussian type with sin- 
gle peak. The continuous shift of the distribution with the temperature over the 
expected critical region is consistent with a second order phase transition in accord 
with the results of the 4 — e expansion method. Final confirmation of the order of 
the transition is done with the values of the critical exponents discussed below. 

3.2 Critical temperature 

Accurate calculation of critical exponents requires a precise determination of the 
inverse critical temperature Kc- An efficient method to determine Kc for a second 



order transition is to measure the Binder cumulant |17| for various system size 
and to locate the cross point in the space of K. On sufficiently large lattices where 
subleading corrections from the finite lattice size L are ignored, the Binder cumulant 
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Ul{K) defined by 



UUK) 



^ 1 (m^) 



3( 



2\2 



1 

X 



becomes independent of L at the transition point Kc [17|: 

and the slope of Ul{K) in K at A'c increases as L becomes large. In Fig.^ are shown 
our results of the Binder cumulant near the crossing point. The values for Ul{K) 
are obtained with the histogram method using the data at K=0.935. 

The deviation from the relation (^) observed in Fig.§ can be explained by the 
finite-size confluent corrections. The leading L'/L-dependence in the deviation of 



the crossing point K* from the critical point Kc is estimated by Binder 1 17] as 

111 

where b = L' /L. 

We plot (1/ln b, 1/K*) for L = 8, 10, 12, and 14 in Fig. |. The errors for l/K* are 
computed from the jackknife errors for Ul{K). The solid lines in Fig. ^ represent the 
results of a linear least-square fit for each L. We find that the correction with 1/ln 6 
is smaller than that of 3d 0(3) Heisenberg model |jl3|. The extrapolation of \/K* to 
the point l/ln6 = for each L gives the values for l/Kc{L): 1/Kc(8) = 1.06841(21), 
1/A"c(10) = 1.06832(26), 1/Kc{l2) = 1.06833(26), l/Kc{U) = 1.06826(34). Ah of 
them are consistent with each other. The mean value of these results is 1/Kc = 
1.06835(13). A similar fit for all L with a common parameter 1/Kc gives the value 
l/Kc = 1.06836(14) which completely agrees with the mean value. We quote here- 
after 

^ = 1.06835(13), (10) 

Kc 

Kc = 0.9360(1). (11) 
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3.3 The critical exponent v 

dUL 



The slope for 



dK 



is known to scale with a critical exponent u as [p^ ] 

K=Kc 



dUL 



dK 

Using the relation 



~ L^l" . (12) 

K=K^ 



we calculate 



^^^-n TJ^\^F^ g ("^'-^) , (^'^) \ rn^ 
-dK-^^- - + T^^T 1 ^ ^ 



at the estimated Kc = 0.9360. In Fig ^ we plot — ^ as a 



dK 

function of L. From the slope of the solid line in this logarithmic plot we find 

- = 1.337(16), (14) 
V = 0.7479(90). (15) 

by a least-square fit. We repeat the analysis by varying Kc within our estimated 
error, 0.9360(1), and find that the results for are completely consistent with the 
result given here. 



The scaling relation (12) requires a sufficiently large L to ignore the subleading 
corrections. In order to test if our values of L is large enough, we repeat the 
fits excluding the data for the smallest size L = 8, and for L = 8 and L = 10. We 
obtain l/u{L = 8 excluded)=l. 344(36) and l/u{L = 8 and 10 excluded)=1.333(51), 
respectively. Because these results are completely consistent with 1/z^ with all data, 
we conclude that L = 8 is sufficiently large to extract scaling properties. 

3.4 The result for - 

V 

The scaling relation of the magnetization (|m|) at K^, is given by 

(|m|)^^ ~ L-PI^ . (16) 

We study the scaling of {\m\)j^^ at Kc = 0.9360 and obtain P/u = 0.5129(7) from 
the slope of the fitted line in Fig. ^. The fits excluding L = 8, and L = 8 and 10 
give the results consistent with this value {(3/v{L = 8 excluded)=0. 5130(15) and 
(5/i'{L = 8 and 10 excluded)=0. 5127(21)). Not like the case of v in the previous 



subsection, we find that the effect of the error of the Kc on the estimate of (3/u is 
larger than the statistical error 0.0007 at Kc = 0.9360: 

-(Kc = 0.9359) --{Kc = 0.9360) = 0.0009 , 

^(Kc = 0.9360) - ^(Kc = 0.9361) = 0.0011. (17) 

Therefore we should use the value 0.0011 for the error of [i/v: 



^ = 0.5129(11). (18) 
Combined with our estimate for i/, we have 

= 0.3836(46) . (19) 

3.5 The result for - 

V 

For K < Kc the susceptibility x is defined by 

X = VK(^m^y (20) 
The scaling relation of x is given by 

{x)k. - L^^" . (21) 

With a similar method as in the previous sections, we obtain for Kc = 0.9360 
7/1^ = 1.9746(15) from the slope of the fitted line in Fig. |6|. Again, the value of j/u 
depends strongly on the choice of Kc- 

-{Kc = 0.9360) -^{Kc = 0.9359) 
^ ^(Kc = 0.9361) - -(Kc = 0.9360) (22) 

V V 

= 0.0038 . (23) 



Therefore we quote 



^ = 1.9746(38) . (24) 



Combined with our estimate of i^, we get 

7 = 1.477(18) . (25) 
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Using our independent results for [3/v and ^/v, we can check the hyperscahng 
relation 

I 1 d 

- + -- -- =0. (26 

We find 

Z./i.s. = 0.0002 ± 0.003 (27) 
that is consistent with zero to O(10~^). 

3.6 Scaling of y*^ and Kyc 

*— ' A.max 

To make a further check of our results for exponents, we study the finite size scaling 
property of the peak of the connected susceptibility x'^'- 

x' = VK[{m')-{\m\f). (28) 

whose maximum value is expected to behave as 

X'^na.-L^^^"^^- (29) 

Here we add a suffix c for the exponent to make clear the way it is defined. Because 
the pseudocritical coupling constant K^c^^^ where g^ts its maximum value is 
found to be slightly off the range of the applicability of the histogram reweighting 
method for the data at (see the discussion in section ^), we carry out new 
simulations at i^' ~ ^Xmax listed in Table |^ determined by a preparatory simulation. 
With the histogram method applied to these new data we estimate accurate values 
for Xmax -ft^Xmaa: (^^^ Table |3|) . From a least-square fit shown in Fig. we obtain 

{-f/u)c = 1.996(8) . (30) 



This value is slightly larger than that from the scaling of x, 1.9746(38), given in (24). 
The same tendency is observed for the 0(3) Heisenberg model [|l^, |T^. Because the 
quality of the fit for x is better than that for Xmax^ quote ( p^ ) for the value of 
7/1/. 

The scaling property of the pseudocritical coupling K^-c^^^ provides us another 
test of our results: 



K-} ~ K^-^ + aL-^/" . (31) 

^ max ^ ' 

10 



Using our estimate 1/v = 1.337, we fit the data with two parameters, ^ and a, 
to obtain 

Kc = 0.9360(2). (32) 
This value is consistent with our Kc from the crossing points of the Binder cumulant. 

3.7 Restriction of the transition point by Q value 

The scaUng relations (|T2D, (p^), and ( p] ) require that the estimated value of is 
close enough to the real transition point. If we fix K^. far from the real transition 
point in these scaling relations the data will not fit them well any more. 
The quality of a least-square fit is determined by the Q value [p^]: 



oo 



Q{x\n)= r'dtl^-Y e-*/2 (33) 



--1 



where is the weighted sum of squared deviations of data from the fit, and n = 
(number of data points) — (number of fit parameters) is the degree of freedom for 
the fit. We may consider that the fitting procedure is appropriate if 0.1 < Q < 0.9. 
If, on the other hand, Q < 0.1 something is wrong: the error of data may be under- 
estimated or the fitting function may be incorrect, and if Q > 0.9 error of data may 
be over-estimated or we have too many fit parameters. 

In the present case, if we fix Kc far from the real transition point, the quality of 
the scaling fits must become low so that the Q- value decreases to a value less than 
0.1. In Fig I the Q-values of our finite-size scaling fits for , {\m\) and x 

plotted as a function of Kc- We find that Q- value for — — ^ is not so sensitive on 

dK 

Kc, while the Q-values for {\m\) and x depend sensitively on Kc- This difference of 

the dependence on Kc between — ^ and (|?ti|), x is the same as that observed for 

dK 

the 0(3) Heisenberg model [|^]. From the condition that Q < 0.1 we have 

0.9359 <Kc< 0.9364. (34) 

This provides us another consistency check of our analyses. The value obtained 
from the crossing point of the Binder cumulant Kc = 0.9360(1) is well included in 
this region. 
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3.8 Comparison with the results of the 4 — e expansion 

The critical exponents obtained in this work are compiled in Table |2| together with 
the values by the 4 — e expansion methodQ. In our results, the exponents 7/1^, (i/v, 
and V are determined independently and a and 6 are calculated using (hyper)scaling 
relations with the value of other exponents. In the results of the 4 — e expansion 
method, r/ and v are estimated independently and other exponents are calculated 
with T] and u. Our results are completely consistent with those of 4 — e with the 
errors reduced to about halves of them. 

4 Conclusion 

We simulated the three-dimensional 0(4) Heisenberg model by applying the single- 
cluster algorithm, which reduces the dynamical exponent to about zero. The his- 
togram reweighting method with high statistics data confirmed that the transition 
is second order for this model. We performed a precise estimation of the critical 
point to get Kc = 0.9360(1) from the crossing point of the Binder cumulant. The 
critical exponents were calculated using finite-size scaling at Kc- The exponents 
obtained, which are summarized in Table H], are completely consistent with those of 
the 4 — e method with the errors reduced to about halves of them. 

We are grateful to Y. Iwasaki and A. Ukawa for valuable discussions and helpful 
suggestions. Numerical calculations were done on HITAC S820/80 at the National 
Laboratory for High Energy Physics (KEK). We would like to thank members of 
KEK for their hospitality and strong support, and members of the Theory Group 
of the Institute of Physics, University of Tsukuba, for innumerable discussions and 
encouragements. This project is in part supported by a Grant-in-Aid of Ministry 
of Education, Science and Culture (No.6206001 and 02402003). 
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L 


simulation point {K) 


sweeps/ 10"^ 








sweeps/ (tc x 10*^) 


8 


0.892 


3000 


31.6 


1.95 


31.6 


95 


Q 
O 


n no cc 
0.935 


oUOU 


CC 1 1 

ol.l 


o /I n 

z.4y 


z4.y 


1 on 


10 


0.935 


4400 


78.7 


2.55 


32.4 


136 


12 


0.910 


4000 


67.3 


1.98 


49.8 


80 


12 


0.935 


6000 


113.2 


2.48 


37.7 


159 


14 


0.912 


1400 


88.0 


1.82 


56.6 


25 


14 


0.935 


6000 


154.5 


2.47 


43.9 


136 


16 


0.920 


3000 


129.5 


1.98 


62.2 


23 


16 


0.935 


5200 


197.4 


2.32 


48.1 


108 


24 


0.926 


1500 


280.9 


1.87 


92.0 


16 


24 


0.935 


2400 


427.1 


2.21 


71.4 


34 


32 


0.928 


1500 


449.5 


1.62 


117.9 


13 


32 


0.935 


5600 


719.1 


1.95 


88.9 


63 



Table 1: Simulation parameters and statistics. (C) is the mean cluster volume; Tc 
is the autocorrelation time in units of sweeps, i.e. in units of the number of cluster 
updates; is the autocorrelation time converted into Metropolis units, i.e. in units 
of updates of whole spins on the lattice: = Tc (C) /V , where V = Li^ \s the lattice 
volume. 
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4-e 


This study 


^ jv = 2 — T] 


1.97(1) 


1.9746(38) 




0.515(5) 


0.5129(11) 


V 


0.73(2) 


0.7479(90) 


7 


1.44(4) 


1.477(18) 


/? 


0.38(1) 


0.3836(46) 


<5 


4.82(5) 


4.851(22) 


a = 2 — dv 


-0.19(6) 


-0.244(27) 



Table 2: Critical exponents of the three-dimensional 0(4) Heisenberg model ob- 
tained by a study with the 4 — e expansion method and by this study. The calcula- 
tion with the 4 — e expansion method is done to seven loops @ and the results are 
quoted in Ref. |||]. In the 4 — e method, independent calculations are done for v and 
T]. In this study, 7/1^, iijv and v are determined independently. Other exponents 
are calculated using (hyper) scaling relations. 
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L 




A.max 


Amax 


8 


0.892 


0.8907(37) 


2.118(09) 


12 


0.910 


0.9109(05) 


4.793(28) 


14 


0.912 


0.9144(06) 


6.577(48) 


16 


0.920 


0.9183(08) 


8.390(46) 


24 


0.926 


0.9253(06) 


18.89(16) 


32 


0.928 


0.9289(01) 


33.97(32) 



Table 3: Results for the pseudocritical coupling Kyc and the maximum values 
Xmax of the connected susceptibility x'^- The K dependence of x'^ is determined 

by the histogram reweighting method using the data simulated at Kgimu on an 
lattice. 
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0.31 0.32 0.33 0.34 0.35 

energy density 

Figure 1: Histogram of the energy density for L=32 near the critical couphng. 
Each distribution is normahzed to unit area. Histograms are for K = 0.929, 0.935, 
and 0.939 from the left to the right, respectively. The histograms for K = 0.929 
and 0.939 are obtained by reweighting the measured histogram at ii'=0.935. 



18 




Figure 2: The Binder cumulant Ul ss & function of the inverse temperature K 
for L = 8, 10, 12, 14, 16, 24, and 32. Ul{K) is computed with the histogram 
reweighting method using the data at K = 0.935. 



19 



1.07 



K 



1.069 - 



1.068 



1.067 - 



1.066 



n — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — rp — r 




012345678 

l/ln(L'/L) 



Figure 3: The crossing point of Binder cumulants Ul{K) and Ul'{K) for L = 
8(squares), lO(triangles), 12 (diamonds), and 14(circles) with different L'. Sohd hnes 
correspond to hnear least-square fits for each L. The critical coupling is estimated 
as Kc = 0.9360(1) by extrapolating these lines to the limit l/ln(L'/L) = 0. 
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Figure 4: Scaling of the slope dUi/dK of the Binder cumulant at Kc = 0.9360(1) 
as a function of the lattice size L. The slope of the solid line given by a linear 
least-square fit leads to an estimate of the critical exponent 1/u = 1.337(16). The 
jackknife errors for dUi/dK are smaller than the size of symbols. 
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Figure 5: Magnetization (|m|) at Kc = 0.9360(1) as a function of the lattice size L. 
A least-square fit gives l3/u = 0.5129(7). The jackknife errors for (|m|) are smaller 
than 1/10 of the size of symbols. 
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Figure 6: Susceptibility x at if c = 0.9360(1) as a function of L. A least-square fit 
gives = 1.9746(15). The jackknife errors for xi^c) are smaller than 1/10 of the 
size of symbols. 
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Figure 7: The maximum height of the connected susceptibility Xmax- ^ least- 
square fit gives {^/i')c = 1.996(8). The jackknife errors for Xmax smaller than 
1/5 of the size of symbols. 
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Figure 8: Q values of least-square fits for dllL/dK (circles), (squares), and 

X (triangles) for various fixed Kc- 
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